The R Spatial Ecosystem

A few main packages are dedicated to spatial data handling like import/export, display, geoprocessing, spatial ststistics.

“Historical” Packages

rgdal: interface between R and GDAL (Geospatial Data Abstraction Library) and PROJ4 libraries.
sp: classes and methods for spatial data in R.
rgeos: interface between GEOS (Geometry Engine - Open Source) library and R: area, perimeter, distances, dissolve, buffer, overlap, union, contains…

These packages are still widely used.

Simple Features for R

sp, rgeos and rgdal functionnalities in one package.
Easier data handling, simpler objects.
Main author and maintainer: Edzer Pebesma (also sp author)
Compatibility with the pipe synthax and tidyverse operators.
First release: October 20, 2016
Website: Simple Features for R

sf objects data structure:

format sf

Using `sf`

Data Import

library(sf)
mtq <- st_read("../data/mtq/martinique.shp")

Reading layer `martinique' from data source `C:\Users\Kim Antunez\Documents\Projets_R\Geovisualisation\flo\data\mtq\martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID):    32620
proj4string:    +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs

Data Display

Default

plot(mtq)

Only geometry

plot(st_geometry(mtq))

Centroids Extraction

mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add = TRUE, cex = 1.2, col = "red", pch = 20)

Distance Matrix

mat <- st_distance(x = mtq_c, y = mtq_c)
mat[1:5, 1:5]

Units: m
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]     0.000 35297.56  3091.501 12131.617 17136.310
[2,] 35297.557     0.00 38332.602 25518.913 18605.249
[3,]  3091.501 38332.60     0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702     0.000  7177.011
[5,] 17136.310 18605.25 20226.198  7177.011     0.000

Polygons Aggregation

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2, border = "red")

Buffer Zone

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col = "lightblue")
plot(st_geometry(mtq_u), add = T, lwd = 2)
plot(st_geometry(mtq_b), add = T, lwd = 2, border = "red")

Polygon Intersection

m <- rbind(c(700015, 1624212), c(700015, 1641586), c(719127, 1641586), c(719127, 
    1624212), c(700015, 1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border = "red", lwd = 2, add = T)

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col = "red", border = "green", add = T)

Voronoi Polygons

google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)

mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join = st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col = "lightblue")

Maps with R

Several solutions are available:

ggplot2 users can have a look to ggplot2 mapping features (geom_sf) that can mix nicely with ggspatial.
For more advanced mapping features in a ggplot2-like syntax have a look to tmap
cartography is based on base graphics and allow most of basic and advanced cartographic representations. Full disclosure: the speaker is the maintainer of cartography.

Here we will focus on cartography and do small examples with ggplot2.

`cartography`

Data Preparation

library(sf)
# Import geo layer
sm <- st_read(dsn = "../data/rhone/seine_maritime.geojson", stringsAsFactors = F)

Reading layer `seine_maritime' from data source `C:\Users\Kim Antunez\Documents\Projets_R\Geovisualisation\flo\data\rhone\seine_maritime.geojson' using driver `GeoJSON'
Simple feature collection with 718 features and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 0.06560943 ymin: 49.25139 xmax: 1.790225 ymax: 50.07121
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

dep <- st_read(dsn = "../data/rhone/dep.geojson", stringsAsFactors = F)

Reading layer `dep' from data source `C:\Users\Kim Antunez\Documents\Projets_R\Geovisualisation\flo\data\rhone\dep.geojson' using driver `GeoJSON'
Simple feature collection with 96 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -5.138996 ymin: 41.36216 xmax: 9.559598 ymax: 51.089
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

# change projection (lambert93)
sm <- st_transform(sm, 2154)
dep <- st_transform(dep, 2154)

# Import dataset
csp <- read.csv("../data/rhone/data_seine_maritime.csv")
# merge geolayer and dataset
sm <- merge(sm, csp, by = "INSEE_COM", all.x = TRUE)

Proportional symbols

library(cartography)

plot(st_geometry(sm))
propSymbolsLayer(sm, var = "act")
title("Active Population")

# Custom map of active population
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
bb <- st_bbox(sm)
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
plot(st_geometry(sm), col = "cornsilk2", border = NA, lwd = 0.5, add = T)
propSymbolsLayer(sm, var = "act", col = "darkblue", inches = 0.5, border = "white", 
    lwd = 0.7, symbols = "square", legend.style = "e", legend.pos = "topleft", 
    legend.title.txt = "Labor Force\n(2014)", legend.values.rnd = 0)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce in Seine-Maritime", sources = "Insee, 2018", 
    author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE, 
    scale = NULL, north = FALSE)

Typology Map

par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
typoLayer(sm, var = "cat", border = "ivory", lwd = 0.5, legend.values.order = c("agr", 
    "art", "cad", "int", "emp", "ouv"), col = c("#e3b4a2", "#a2d5d6", "#debbd4", 
    "#b5dab6", "#afc2e3", "#e9e2c1"), add = TRUE, legend.pos = "n")

legendTypo(title.txt = "Dominant Category", col = c("#e3b4a2", "#a2d5d6", "#debbd4", 
    "#b5dab6", "#afc2e3", "#e9e2c1"), categ = c("Agriculteurs", "Artisans", 
    "Cadres", "Prof. Inter.", "Employés", "Ouvriers"), nodata = F)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce Distribution in Seine-Maritime", sources = "Insee, 2018", 
    author = "Kim & Tim, 2018", col = "darkblue", coltitle = "white", frame = TRUE, 
    scale = NULL, north = FALSE)

Choropleth Map

sm$pcad <- 100 * sm$cad/sm$act

bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
choroLayer(sm, var = "pcad", breaks = bks, col = cols, border = "grey80", lwd = 0.4, 
    legend.pos = "topleft", legend.title.txt = "Share of managers\n(%)", add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

library(cartogram)
sm_c1 <- cartogram_cont(sm, "act", 3)
sm$pcad <- 100 * sm$cad/sm$act

bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)
cols <- carto.pal("green.pal", 3, "wine.pal", 3)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
choroLayer(sm_c1, var = "pcad", breaks = bks, col = cols, border = "grey80", 
    lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Share of managers\n(en %)", 
    add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

Gridded Map

These maps may allow to hide or, at least, diminish the MAUP.

sqrt(quantile(st_area(sm), 0.9))

3904.689 m

grid <- getGridLayer(x = sm, cellsize = 3900 * 3900, type = "regular", var = c("cad", 
    "act"))
grid$pcad <- 100 * grid$cad/grid$act
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
choroLayer(grid, var = "pcad", breaks = bks, col = cols, border = "grey80", 
    lwd = 0.4, legend.pos = "topleft", legend.title.txt = "Share of managers\n(en %)", 
    add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

Smoothed Map

grid$cad100 <- grid$cad * 100
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(dep), col = "ivory", border = "ivory3", xlim = bb[c(1, 3)], 
    ylim = bb[c(2, 4)])
smoothLayer(x = grid, var = "cad100", var2 = "act", typefct = "exponential", 
    span = 4000, beta = 2, breaks = bks, col = cols, legend.pos = "topleft", 
    mask = st_buffer(sm, 0), legend.title.txt = "Share of managers\n(%)*", border = "grey90", 
    lwd = 0.2, add = T)
layoutLayer(title = "Managers", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
    theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center", 
    frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)
text(x = 518000, y = 6921123, font = 3, cex = 0.8, labels = "*Lissage par potentiels\n fonction exponnentielle\n span = 2.5km, beta = 2")

`ggplot2`

Proportional symbols layer:

library(ggplot2)

ggplot() +
  geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
  geom_sf(data = sm %>%  st_centroid(),
          
          aes(size= cad), colour="#E84923CC", show.legend = 'point') +
  scale_size(name = "Active population",
             breaks = c(100,1000,5000),
             range = c(0,20)) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(sm)[c(1,3)],
           ylim = st_bbox(sm)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(title = "Active population",
    caption = "Insee, 2018\nKim & Tim, 2018")

Choropleth layer:

ggplot() +
  geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
  geom_sf(data = sm, aes(fill = pcad), colour = "grey80") +
  scale_fill_gradientn(name = "Share of managers\n(%)",
                       
                       colours = carto.pal("green.pal", 3,"wine.pal",3),
                       values=bks/max(bks))+
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(sm)[c(1,3)],
           ylim = st_bbox(sm)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(title = "Managers",
       caption = "Insee, 2018\nKim & Tim, 2018")

Exploration of sirene & osm db for restaurants

CRAN Task View: Analysis of Spatial Data

See file data_prep.R for data extraction.

# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)

Display Points

# load data
adm <- readRDS("../data/iris_p13.rds")
feat_sir <- readRDS("../data/sir_p13.rds")
feat_osm <- readRDS("../data/osm_p13.rds")

nrow(feat_sir)

[1] 1066

nrow(feat_osm)

[1] 585

par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))

plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, frame = FALSE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, cex = 0.7, 
    col = "red")

plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add = T, cex = 0.5)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T, 
    frame = FALSE)

Count Points in Units

In admin units

inter_osm <- st_intersects(adm, feat_osm)
inter_sir <- st_intersects(adm, feat_sir)

adm <- st_sf(adm[, 1:2, drop = T], n_osm = sapply(X = inter_osm, FUN = length), 
    n_sir = sapply(X = inter_sir, FUN = length), geometry = st_geometry(adm))

par(mar = c(0, 0, 1.2, 0), mfrow = c(1, 2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_sir", inches = 0.1)
layoutLayer(title = "SIR", sources = "", author = "", scale = NULL, frame = FALSE)

plot(st_geometry(adm), col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_osm", inches = 0.1)
layoutLayer(title = "OSM", sources = "", author = "", scale = 0.5, north = T, 
    frame = FALSE)

# plot(adm$n_sir, adm$n_osm, asp = T, xlim = c(0,55), pch = 21, bg = 'red')
# text(adm$n_sir, adm$n_osm, labels = adm$CODE_IRIS, cex = 0.5, pos = 2) x
# <-lm(formula = adm$n_osm~adm$n_sir) abline(a = x$coefficients )

In a regular grid

fufun2 <- function(feat, adm, cs = 500) {
    grid <- st_make_grid(x = adm, cellsize = cs, what = "polygons")
    x <- st_intersects(grid, feat)
    gg <- st_sf(n = sapply(X = x, FUN = length), grid)
    plot(adm$geometry, col = "ivory1", border = "ivory3", lwd = 0.5, bg = "#FBEDDA")
    choroLayer(gg[gg$n > 0, ], var = "n", border = NA, breaks = c(1, 4, 6, 9, 
        12, 14, 17, 19, 22), col = carto.pal("wine.pal", 8), add = TRUE)
}

par(mar = c(0, 0, 1.2, 0), mfrow = c(2, 3))
fufun2(feat_sir, adm, 50)
fufun2(feat_sir, adm, 100)
fufun2(feat_sir, adm, 150)
fufun2(feat_osm, adm, 50)
fufun2(feat_osm, adm, 100)
fufun2(feat_osm, adm, 150)

Interactive Visualisation

Explore datasets interactively

The default is not bad:

mapview(feat_sir) + mapview(feat_osm, col.regions = "red")

mapview allow hightly customizable maps:

mapview(feat_sir, map.types = "OpenStreetMap", col.regions = "#940000", label = paste(feat_sir$L2_NORMALISEE, 
    feat_sir$NOMEN_LONG, sep = " - "), color = "white", legend = TRUE, layer.name = "SIRENE", 
    homebutton = FALSE, lwd = 0.5) + mapview(feat_osm, col.regions = "#000094", 
    color = "white", legend = TRUE, label = feat_osm$name, lwd = 0.5, layer.name = "OSM", 
    homebutton = FALSE)

These maps, as appealing as they seem to be, are not really suitable for displaying geostatistical info. For data exploration they are really useful.

Smooth Density Analysis

Define a function that compute smoothed restaurant density and map it. KDE method

fufun <- function(feat, adm, title, sigma = 100, res = 50) {
    bb <- as(feat, "Spatial")
    bbowin <- as.owin(as(adm, "Spatial"))
    pts <- coordinates(bb)
    p <- ppp(pts[, 1], pts[, 2], window = bbowin)
    ds <- density.ppp(p, sigma = sigma, eps = res)
    rasdens <- raster(ds) * 1000 * 1000
    proj4string(rasdens) <- "+init=epsg:2154"
    rasdens <- rasdens + 1
    bks <- getBreaks(values(rasdens), nclass = 12, method = "equal")
    cols <- (mapview::mapviewGetOption("raster.palette"))(length(bks) - 1)
    plot(adm$geometry, col = NA, border = NA, main = "", bg = "#FBEDDA")
    plot(rasdens, breaks = bks, col = cols, add = T, legend = F)
    legendChoro(pos = "topright", cex = 0.7, title.cex = 0.7, title.txt = "resto per km2)", 
        breaks = bks - 1, nodata = FALSE, values.rnd = 0, col = cols)
    layoutLayer(title = paste0(title), scale = 0.5, tabtitle = TRUE, frame = FALSE, 
        author = "Map data © OpenStreetMap contributors, under CC BY SA. SIRENE 2018", 
        sources = "cartography 2.0.2")
    return(invisible(rasdens))
}

par(mar = c(0, 0, 1.2, 0), mfrow = c(2, 3))
fufun(feat_sir, adm, "SIR 50", 50, 20)
fufun(feat_sir, adm, "SIR 100", 100, 20)
fufun(feat_sir, adm, "SIR 250", 250, 20)
fufun(feat_osm, adm, "OSM 50", 50, 20)
fufun(feat_osm, adm, "OSM 100", 100, 20)
fufun(feat_osm, adm, "OSM 250", 250, 20)

reproducibility

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mapview_2.4.0       maptools_0.9-2      raster_2.6-7       
 [4] sp_1.3-1            spatstat_1.56-1     rpart_4.1-13       
 [7] nlme_3.1-137        spatstat.data_1.3-1 ggplot2_3.0.0      
[10] cartography_2.1.1   sf_0.6-3            rmdformats_0.3.3   
[13] knitr_1.20         

loaded via a namespace (and not attached):
 [1] viridisLite_0.3.0    foreach_1.4.4        R.utils_2.6.0       
 [4] shiny_1.1.0          assertthat_0.2.0     stats4_3.5.1        
 [7] highr_0.7            yaml_2.2.0           pillar_1.2.3        
[10] backports_1.1.2      lattice_0.20-35      glue_1.3.0          
[13] digest_0.6.15        promises_1.0.1       polyclip_1.9-1      
[16] colorspace_1.3-2     htmltools_0.3.6      httpuv_1.4.4.1      
[19] Matrix_1.2-14        R.oo_1.22.0          plyr_1.8.4          
[22] pkgconfig_2.0.1      questionr_0.6.3      bookdown_0.7        
[25] purrr_0.2.5          xtable_1.8-2         webshot_0.5.0       
[28] scales_0.5.0         tensor_1.5           satellite_1.0.1     
[31] later_0.7.3          spatstat.utils_1.9-0 tibble_1.4.2        
[34] mgcv_1.8-24          withr_2.1.2          lazyeval_0.2.1      
[37] gdalUtils_2.0.1.14   magrittr_1.5         mime_0.5            
[40] deldir_0.1-15        evaluate_0.11        R.methodsS3_1.7.1   
[43] foreign_0.8-70       class_7.3-14         tools_3.5.1         
[46] formatR_1.5          stringr_1.3.1        munsell_0.5.0       
[49] bindrcpp_0.2.2       compiler_3.5.1       e1071_1.6-8         
[52] rlang_0.2.1          classInt_0.2-3       units_0.6-0         
[55] grid_3.5.1           iterators_1.0.10     rstudioapi_0.7      
[58] htmlwidgets_1.2      goftest_1.1-1        crosstalk_1.0.0     
[61] miniUI_0.1.1.1       base64enc_0.1-3      labeling_0.3        
[64] rmarkdown_1.10       gtable_0.2.0         codetools_0.2-15    
[67] abind_1.4-5          DBI_1.0.0            R6_2.2.2            
[70] dplyr_0.7.6          rgdal_1.3-3          rgeos_0.3-28        
[73] bindr_0.1.1          rprojroot_1.3-2      stringi_1.1.7       
 [ reached getOption("max.print") -- omitted 6 entries ]

Geovisualization and Spatial Analysis